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Abstract. 

A detailed analysis of dynamics of cosmological models based on R n gravity is 
presented. We show that the cosmological equations can be written as a first order 
autonomous system and analyzed using the standard techniques of dynamical system 
theory. In absence of perfect fluid matter, we find exact solutions whose behavior 
and stability are analyzed in terms of the values of the parameter n. When matter 
is introduced, the nature of the (non-minimal) coupling between matter and higher 
order gravity induces restrictions on the allowed values of n. Selecting such intervals 
of values and following the same procedure used in the vacuum case, we present exact 
solutions and analyze their stability for a generic value of the parameter n. From 
this analysis emerges the result that for a large set of initial conditions an accelerated 
expansion is an attractor for the evolution of the R n cosmology. When matter is 
present a transient almost-Friedman phase can also be present before the transition to 
an accelerated expansion. 
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1. Introduction 

At the present time, there are many good reasons to consider General Relativity (GR) 
as our best theory for the gravitational interaction. It accounts for the shortcomings 
of Newtonian gravity and fits beautifully all the tests that we have been able to devise 
at laboratory and Solar System levels. Nevertheless, in the last few decades, there is 
more and more evidence that GR may be incomplete. In particular, the difficulties 
to match GR with quantum theory, the flattening of galactic rotation curves and, 
more recently, the cosmic acceleration can be considered hints that some aspects of 
gravitational interaction are still unknown. 

It is for these reasons that much recent effort has been spent on the study of 
generalizations of the Einstein theory and in particular the so-called Extended Theory 
of Gravitation (ETG). 

One of the most interesting classes of ETG, called non-linear gravity theories or 
higher-order theories of gravity, is based on gravitational actions which are non-linear 
in the Ricci curvature (R) and/or contain terms involving combinations of derivatives 
of R p] EI] . The higher order theories of gravity arise in a wide range of different 
frameworks. For example, it can be shown that for quantum field theory in curved 
spacetime, the renormalization of the stress energy tensor implies the introduction of 
higher order corrections to the GR gravitational action In addition, the low energy 
limit of D = 10 superstring theory is a theory of gravity with additive higher order 
corrections to the Hilbert-Einstein (HE) term E] • Finally, the general form of the 
vacuum action for Grand Unification Theories (GUT) is a higher order theory and the 
GUT models obtained by relaxing the constraint of linear gravity have, in general, better 
physical properties than the standard GR [Zj. 

Unfortunately most of the efforts to achieve an insight into the physics of these 
models has been frustrated by the complexity of the field equations. In fact, the 
presence of non-linearities and higher order terms makes it hard to achieve both exact 
and numerical solutions which can be compared with observations. Such difficulties, 
added to other problems connected with the presence of non-gaugeable ghosts in the 
particle spectrum of these theories jH] , have led researchers to pay more attention towards 
other alternative gravity theories such as braneworld models [H]. 

The recent realization that the universe is currently undergoing an accelerated 
expansion phase and the quest for the nature of the quintessence field has renewed the 
interest in higher order theories of gravity and their relation to cosmology. This is due to 
the fact that the higher order corrections to the HE action can be viewed as an effective 
fluid and this fluid can emulate the action of the homogeneous part of the quintessence 
field PU]- Hence, in this Curvature Quintessence scenario, what we observe as a new 
component of the cosmic energy density is an effect of higher order corrections to the 
HE Lagrangian. This approach has several advantages over the standard quintessence 
scenario. For example, we do not need to search for a quintessence scalar field and such 
theories are more easily constrained by observations. 
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A key element of the curvature quintessence scenario is the choice of the form 
of the Lagrangian density C for the gravitational interaction. A simple choice is a 
generic power law of the Ricci scalar R. This R" -gravity model has been shown to be a 
good candidate for describing the present evolution of the Universe [TTH ITU FEU ITIj. 
Comparison with data coming from supernovae data and the age of the universe as 
measured by WMAP [15 j shows that, at least in the vacuum case, these models are 
both in agreement with the observations and allow accelerating solutions ^H]- This 
result is very important, but incomplete. In fact, we do not know what the stability of 
these solutions are and which role they play in the global behavior of the underlying 
cosmological model. In addition, when perfect fluid matter is introduced, the equations 
become even more complex, making it extremely difficult to find exact solutions. 

These problems can be (at least partially) addressed using the dynamical system 
approach. This powerful method was first developed by Collins and then by Ellis 
and Wainwright (see J7j for a wide class of cosmological models in the GR context). 
Some work has also been done in the case of scalar fields in cosmology and for scalar- 
tensor theories of gravity [TBI 113 1211] • Studying cosmologies using the dynamical 
systems approach has the advantage of offering a relatively simple method to obtain 
exact solutions (even if these only represent the asymptotic behavior) and to obtain a 
(qualitative) description of the global dynamics of these models. Such results are very 
difficult to obtain by other methods. 

The purpose of this paper is to apply the dynamical system approach to higher order 
gravity and in particular to the R n — gravity models. The paper is organized as follows: 
In section |2J we present the basic equations of the model. In section |3J we analyze 
dynamics of R n — gravity in the vacuum case, find exact solutions and determine their 
stability. Section 0] gives a generalization of these results to the matter case. Finally, 
section |S] contains a summary of the results and conclusions. 

Unless otherwise specified, we will use natural units (h = c = ks = 87rG = 1) 
throughout the paper, Greek indices run from to 3 and Latin indices run from 1 to 
3. The symbol V represents the usual covariant derivative and d corresponds to partial 
differentiation. We use the (+, — , — , — ) signature and the Riemann tensor is defined by 

rp, =w p x - w p , + w a we + w^w p „ (i) 

where the is the Christoffel symbol (i.e. symmetric in the lower indices), defined 

by 

The Ricci tensor is obtained by contracting the first and the third indices 

R^u = g pX R Pf i\u ■ (3) 
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2. Basic Equations 

If we relax the assumption of linearity in the gravitational action, the general coordinate 
invariance allows infinitely many additive terms to the HE action: 

S = J d 4 x^ {A + cqR + cii? 2 + c 2 R tll/ R> xu . . . + C M ) , (4) 

where Cm is the matter Lagrangian and we have included corrections up to the fourth 
order (the fourth order term R llua pR llva P has been neglected as a consequence of the 
Gauss-Bonnet theorem. 

The action (J3J) is not canonical because the Lagrangian function contains derivatives 
of the canonical variables of order higher than one. This means that, not only do we 
expect higher order field equations, but also the validity of the Euler-Lagrange equations 
is compromised. 

This problem is particularly difficult in the general case, but can be solved for 
specific metrics. In fact, once the components of the metric tensor are specified, we can 
define a new set of generalized coordinates in such a way to reduce formally the degree 
of the Lagrangian. In these new variables, the higher order field equations correspond 
to a system of second order differential equations together with a constraint equation 
associated with the definition of the variables themselves [22]- As will be clear in the 
following, this procedure is crucial in order to both write the field equations for R n - 
gravity and to apply the dynamical system approach. 

In homogeneous and isotropic spacetimes, the Lagrangian in can be further 
simplified. Specifically, the variation of the term R^R^ can always be rewritten in 
terms of the variation of R 2 Thus, the "effective" fourth order Lagrangian in 

cosmology contains only powers of R and we can suppose, without loss of generality, 
that the general form for a non-linear Lagrangian will be a generic function of the scalar 
curvature: 

L = ^g-(f(R) + C M ) . (5) 
By varying $ equation (jSJ), we obtain the fourth order field equations 

f'(R)R, u - \f{R)g»v = f'(R)' aP (g a »gp„ - g^) + f£ , (6) 



where T = and the prime denotes a derivative with respect to R. It 

v-g og^ v 

is easy to check that standard Einstein equations are immediately recovered if f[R) = R- 
When f'(R) 7^ the equation (jHJ) can be recast in the more expressive form 

1 

T 



Gfii, — R^ v ^g^ v R — — T^ u + T^ u , (7) 



where 



t% = jj^ if(R) - Rf'(R)] + f'(R)^(g«^ - g a pg^)\ 



\ In the present work we will adopt a "metric approach to derive the field equations. Alternatively, one 
can use the "Palatini approach" in which the fields "<?" and "W" are considered independent [111 I14| . 
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and 

rpM _ 1 rpM fq\ 

is an effective stress-energy tensor for standard matter. This step is conceptually very 
important: we have passed from a model in which gravity has a complicated structure 
to a model in which the gravitational field has the standard GR form with a source 
made up of two fluids: perfect fluid matter and an effective fluid (curvature fluid) that 
represents the non-Einsteinian part of the gravitational interaction §. 

If we adopt the above picture, there is an important issue related to the conservation 
properties of the two new fluids. In fact, looking at the Bianchi identities for equation 
(J7J), we expect a conservation equation with interaction terms between the curvature 
fluid and standard matter. However this is not the case. In fact, writing down the 
Bianchi identities we have 



f) _ rpTOT,U rpR-u _ f (R) T3,VrpM , \ rpM]U /in\ 

■ L ii,v L iiv f'^(R) v f'(R) ^ ' 



Using the field equations and the definition of the Riemann tensor, it is easy to show 
that the sum of the first two terms of the RHS of the previous expression is zero. Thus, 
T™ T ' ,U oc T^l and the total conservation equation reduces to the one for matter only. A 
general, independent proof of this result has been given by Eddington |21| and then by 
others (like in [25]). They showed that the first variation for the gravitational action is 
divergence free regardless of the form of the invariant that we choose for the Lagrangian. 
This means that no matter how complicate the effective stress energy tensor T™ T is, 
it will always be divergence free if T^ v = 0. As a consequence, the total conservation 
equation reduces to the one for standard matter only. This is a very important feature 
of the curvature fluid and it is the fundamental reason why we can use the Friedmann 
metric in these models even though the curvature fluid may in general possess heat flux 
and anisotropic pressure 

This paper will focus on a particular higher order model: _R ra -gravity, whose action 
for the gravitational interaction reads 

A = j d*xV=g[x(n)R n + L M } , (11) 

where x{ n ) is a function of n that reduces to 1 for n = 1. Whatever the value of n 
we will suppose that the sign of x is k e Pt positive in such a way that the attractive 
character of the gravitational interaction is not compromised. 
The field equations for this theory read 

nR n ~ l G, v = x(n)- 1 ^ + ^(1 - n)R n + [n{n - l)R n ~ 2 R^ 
+n ( n - l)( n - 2)R n - 3 R a R?] (g^ gf3u - g aP g^) (12) 

§ Obviously we do not need to restrict ourselves to the f(R) Lagrangians to make this transformation. 
In general, we could have started with (@J at whatever order and define an effective stress energy tensor 
in such a way that the field equations have the form 10. 
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and for R ^ they can be recast as 



rpM | 
iiu J- 



fn- 1' 



/ -,\/ „.R ,a R^ 
+(n - l)(n - 2) 



7? 

(SVS^ - 9ap9w) ■ (13) 



i? 2 

In this way, the non-Einsteinian part of the gravitational interaction can be modelled 
as an effective fluid which, in general, posses thermodynamic properties different from 
standard matter. 

In the Friedmann-Lemaitre- Robertson- Walker (FLRW) metric, the above equations 
can be written as 

2n- + n(n - 1)H— + n(n - 1)— + n(n -l)(n- 2)^- - (1-n) — 
a v J R K J R v A 1 R 2 v 7 3 

q + 3w) = , (14) 



with 



3nx{n)R n - 1 

H 2 + ^- + H-(n-l)-—(l-n) r = 0. (15) 

a 2 R K ' 6n y ' 3nx(n)R n ~ l v ; 

(a a 2 K \ , . 

R=-6[~ + — + — ) 16 



where H = a/a, k is the spatial index and we have assumed the standard matter to be 
a perfect fluid with a barotropic index w. Note that in the above equations we have 
considered a and R as two independent fields. Such a position is standard in canonical 
quantization of higher order gravitational theories |2Z] and is very useful in this case 
because it allows to write the equations (fT2"|) as a system of second order differential 
equations. 

The system (jl4H16|) is completed by the conservation equation for matter, which as 
we already pointed out, is the same as standard GR, i.e. 

p + 3Hp(l + w) = . (17) 

We start first with an analysis of the vacuum case and then present the general case in 
which matter is present. 

3. The Vacuum Case 

Equations (jl4H15|) are the starting point of a dynamical analysis of cosmological models 
based on R n gravity. In vacuum (p = 0), their form suggests the following choice of 
expansion normalized variables: 

* = ^("-i)> 
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The variable x is associated with the rate of variation of the Ricci curvature, y represents 
a measure of the expansion normalized Ricci curvature and K the spatial curvature of 
the Friedmann model. Following [2B|, if we consider the Ricci curvature as a scalar field, 
we can think of x as the kinetic term for this field and y as its potential. However, this 
analogy does not work completely, because x appears only linearly in the Friedmann 
equation. 

Defining the new time variable M = In a, the cosmological equations are equivalent 
to the autonomous system 

x' = 2 + x - x 2 + (2 + nx) + (2 + x)K , 

1 — n 

v > = ^JL + 2v (J^ + 2 + k ) , (19) 

n — 1 \\ — n J 

K > = 2K{1 + + K) , 

1 — n 

where the prime represents the derivative with respect to our new "time" coordinate 
and the dynamical variables are constrained by 

l + K + x-y = 0. (20) 

It is important to stress that even if the cosmological equations and the dynamical 
variables are defined for every value of n the dynamical equations (fTT?|) are only valid 
for n 7^ 1. This is due to the procedure we have used to obtain this system. However, 
since the case n = 1 corresponds to standard GR, we will limit our discussion ton^ 1. 

3.1. Finite analysis 

The system (II 9|) can be further simplified because the constraint ()20|) allows to write 
the equation for K as a combination of the other two variables x and y. In this way, we 
can reduce to the two dimensional system 

i (2 + nx)y , , 

x' = (2 + x )y-2x-2x 2 -^— J — , 21 

n — 1 

V = ^37 [(3 - 2n)x -2y + 2(n - 1)] , 

and use (J2(J|) to determine the value of K. Note that if y — the above system implies 
y' = and the x axis is an invariant submanifold for the phase space. This means that 
if the initial condition for the cosmological model is y ^ a general orbit can approach 
y = only asymptotically. As a consequence, there is no orbit that crosses the x axis 
and no global attractor can exist. 

Setting x' = and y' = 0, we obtain five fixed points 

K = -l, 
K = , 

2(71-2)1 

irrr}- K = °- < 22 » 

+ -2(n-l)}, K = 2n 2 - 2n - 1 , 



A = 


{y- 


-> 0,2 - 




B = 


{y- 


-> 0, x - 




C = 




An - 


5 


{,. 


2n- 








V = 


{y- 


-> 2{n- 


■l) 2 ,x 
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Table 1. Coordinates of the fixed points, eigenvalues, and solutions for i? n -gravity in 



Point Coordinates (x,y) 



Eigenvalues 



A 
B 
C 

V 



[0,0] 
[-1,0] 

2(n-2) 4n-5 
2n-l ' 2n-l 



[2(1 -n) ,2(n-i; 



[-2,2] 

[9 4n-5 1 
L Z ' n-1 J 
5— 4n 4n+2-4ri 2 
n-1 ' 2n 2 -3n+l 

n-2 + ^3n(3n - 4) , n - 2 - ^3^(3™ - 4) 



Point 



Solution 



A 
B 



V 



a = aot 

a = a {t- t ) 1/2 (only for n = 3/2) 

(l-n)(2n-l) 

a = a t «-2 



a 



- if MO 



2n 2 -2n- 

a = a t if = 



with A being a double solution. We can distinguish two classes of fixed point. In 
the case of the first one (A,B), the coordinates are independent of the value of n and 
common to all the R n theories. In the second one (C,T>), the coordinates of the point 
are n-dependent and their position changes with the parameter n. In particular, the 
coordinates of the point C diverge for n = 1/2. 

Merging occurs for the points C and B at n = 5/4 and for C and T> at n = \{1 + v^3) 
and n = |(1 — y/3). The points A, D would merge at n = 1, but we will not consider 
this case because our system is not defined for this value of n. 

The stability of the fixed points can be obtained by a local stability analysis [3Ti] 
and this is summarized in Table |21 above. 

In A the linearized matrix has a positive and a negative eigenvalue which is 
independent of the value of n i.e. it is always a saddle. 

The point B has instead only one fixed (positive) eigenvalue and a second eigenvalue 
which depends on n. Consequently, it is an unstable node for n < 1 and n > 5/4 and a 
saddle otherwise. 

In the case of C, both the eigenvalues depend on n. Yom < — >/3), 1/2 < n< 1 
and n > |(1 + \/o) it corresponds to stable node and for 1 < n < 5/4 is represents an 
unstable node. For the other values of n, it is a saddle point. 

Similarly for V, both the eigenvalues depend on n. For n < V3), n > |(1+V3) 
it is a saddle point, for |(1 — a/3) < n < and 4/3 < n < |(1 — \/3) it is a stable node 
and for the other values of n the eigenvalues are complex so T> represents an attractive 
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Table 2. Stability of the fixed points for R" -gravity in vacuum. We refer to the 
attractive foci as "spirals " 



Point 


n < |(1 - V3) 


|(1 - V3) < n < 


< n < 1/2 


1/2 < n < 1 


A 


saddle 


saddle 


saddle 


saddle 


B 


repulsive 


repulsive 


repulsive 


repulsive 


C 


attractive 


saddle 


saddle 


attractive 


V 


saddle 


attractive 


spiral 


spiral 



Point 


1 < n < 5/4 


5/4 < n < 4/3 


4/3 <n< |(1 + V3) 


n > |(1 + v^) 


A 


saddle 


saddle 


saddle 


saddle 


B 


saddle 


repulsive 


repulsive 


repulsive 


C 


repulsive 


saddle 


saddle 


attractive 


V 


spiral 


spiral 


attractive 


saddle 



focus||. 

A pictorial representation of the phase space in the relevant intervals for n is given in 
figures (1-8). It is clear that in this plane the x axis represents an invariant submanifold 
because no orbit crosses it. 

The coordinates of the fixed points ()22j) can be used to find exact solutions for this 
cosmological model. In fact, by substituting the definitions f|T8|) in the field equations 
we obtain 

H = (x c + - 1 )h 2 , (23) 



n — 1 

which links the coordinates of the fixed point to the form of H(t). If n ^ 1 and 

(n-l)(x c -l)+yc^0 (24) 
this equation can be integrated to give 

(Vc \ 
I — xc ■ (25) 
n-lj 

For the point A we have 

a = a (t - t ) , (26) 

which represent a Milne evolution and is independent of the value of the parameter n. 
For the point B we have 

a = a (i-i ) 1/2 , (27) 

|| The values of n for which the stability of the fixed points changes are called bifurcations for the 
system (j2U- In what follows we will not deal with these specific cases leaving this discussion for a 
future work. 
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but the direct substitution in the equation reveals that this solution is valid only for 
n = 3/2. For the point C we have 

(l-n)(2«-l) 

a = a {t- 1 ) - 2 , (28) 

which corresponds to a solution obtained elsewhere ^U] by direct integration of the field 
equations. The difference is that the dynamical system approach used here allows us to 
deduce the stability of this solution. This issue is particularly important when we use 
the constraint on n from observations ^E]. For particular values of this parameter, C 
is not only an attractor but also the only finite attractor for the phase space. For the 
point T> we have again a Milne evolution: 

a = a (t- t ) , (29) 

but this time the value of n influences the behavior of the solution. In fact for n ^ 1=t: 2 v/ ^ 
the constant Oq is given by 

a °- 2„3-2„.-l ' < 30 > 
otherwise the constant a remains undetermined. It follows that the value of n 
determines the evolution rate unless n = 1:b 2 ■ 

It might seem surprising that these solutions posses less than the four integration 
constants which we would expect when solving a fourth order differential equation. This 
fact can be explained by noting that the solutions for the fixed points are asymptotic 
in nature so we can always imagine that the other constants occur in a transient part 
of the solution that is not relevant in the asymptotic regime. 

It is also useful to define the deceleration parameter q in terms of the dynamical 
variables. We have 

q = —=2 - 1 = ~ x 7 31 

H z n — 1 

which means, for q > 0, 

71 > 1 { n < 1 (32) 

x(n-l)+y<0 1 x(n-l)+y>0 

Thus the phase space is divided into two regions which correspond, once n is fixed, 
to accelerating and decelerating evolution. It is clear that these conditions are never 
satisfied for the points A, £>, T>. For the point C the constraint on n are 

n > |(1 + a/3) expansion, 

| < n < 1 expansion, (33) 
n < \{1 — a/3) contraction. 

3.2. Asymptotic analysis 

Since the phase plane is not compact it is possible that the dynamical system (|21j) 
has a nontrivial asymptotic structure. Thus the above discussion would be incomplete 
without checking the existence of fixed points at infinity (i.e. when the variables 
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x or y diverge) and calculating their stability. Physically, such points represents 
regimes in which one or more of the terms in the Friedmann equation (fT5|) become 
dominant. For example, if x — > oo, the "kinetic" part of the curvature energy density 
dominates over the "potential" part. The asymptotic analysis can easily be performed 
compactifying the phase space using the so called Poincare method. In the vacuum case, 
the compactification can be achieved by transforming to polar coordinates (r(AT), 9(M)): 

x — » rcos(#) , y — » rsin(#) , (34) 

and substituting r — > so that the regime r — ► oo corresponds to 1Z — > 1. Using the 
coordinates (jMj) and taking the limit 1Z — > 1, the equations (|21|) can be written as 

(9 - 8n) cos(fl) - cos(3fl) - 7 sin(0) + sin(3#) 

1^ ~~ * 77 77 (.35J 

4(n — 1) 



cos(fl) sin(0)[sin(e) -cos(„,j 

(n-l)(l-7i) ' (36) 

Starting from (|35H36j) we can analyze the asymptotic features of the system (J2Tj) in the 
same way as we did in the finite case. Indeed, since equation (|33j) does not depend on 
the coordinate 71 we can find the fixed points using only equation (jHEJ). Setting 9' = 0, 
we obtain six fixed points which are given in Table ©■ Their angular coordinate does 
not depend on the value of the parameter n, so they are the same irrespective of the 
power of R appearing in the initial Lagrangian. 

The solutions associated with these fixed points are summarized in Table El The 
points Axm^oo and Coc-T'oo are characterized respectively by the kinetic part or the 
potential part of the energy density of the curvature fluid dominating the dynamics. 
Let us derive the solution for Aoo- This point corresponds to x — >• ±oo, y — > 0. In the 
limit x — >• oo the first equation of the system (J2Tj) and (J23*)) reduce to 

x -> -2x 2 (37) 

and 

— = x . (38) 
Integrating (|3T|) and substituting in ()38|) we obtain 

2 



(39) 



which represents a scale factor that approaches a constant when £ — > £ (-A/" — » A/^). 
The same procedure can be applied to the points P^, C^, JF^ obtaining the same result. 

The points and JF^ represent a state in which the total energy density of the 
curvature fluid dominates and corresponds to closed Einstein universes. For these points 
x and y goes to infinity at the same rate (9 = 7r/4). In this limit, the system (fT9"j) reduces 
to the equation 

1 — 2n 9 . . 

a/ = x 2 40 

n — 1 
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Table 3. Asymptotic fixed points, 9 coordinates and solutions for i? ra -gravity in 



vacuum. 


Point 


e 






Behavior 











w- 


Noo\ = [Ci±f (t- 


tof 


Boo 

Coo 

Voo 


vr/4 
tt/2 

7T 


W 


-Moo 

w- 

w- 


- [Ci ± C 2n _ 1 (t 

JVQ = [d±f (t- 

AU = [Ci±f (*- 


2n-l 

-*o)] 

to)]" 

to)] 2 


c 'oo 

T 

J oo 


5tt/4 
3tt/2 


W 


-Moo 

w- 


- [Ci ± C 2n _J (t 
AToo| = [Ci±f (t- 


2n-l 

-t )] - 1 
to)] 2 



which admits the solution 



n — 1 



(2n-l)(JV-JV 0O ) ' 

where J\f, is an integration constant. In the same way, the (J23)) can be written as 

H' nx 

~H ~ n - 1 ' 
Substituting (}4~Tj) in (|4^j) and integrating we obtain 

\H\ = \H\ = Co\M-Moo\^ , 



(41) 



(42) 



(43) 



where C is a positive integration constant. Integrating again with respect to the variable 
t we obtain 



Ci±Co 



n — 1 



/; 



(44) 



2n - 1 

This solution represents an expansion if n < 1/2 and n > 1, a contraction ifl/2<n<l. 

The value of n influences also the stability of these fixed points (see Table EJ. In 
particular, we see that only the points Boo, T^oo, Foo can be attractors. For the values of 
n suggested by the observations only the points T>oo or Too are attractors and we can 
conclude by saying that, for these values of n, the evolution is always given by (jSJ). 

Summarizing the above results, we can say that the global behavior of vacuum 
R n gravity exhibits transient Milne and power law states that eventually approach the 
solutions C, V or B^, P^, depending on the value of the parameter n and the 
initial conditions. For n far from one, only the points C and T> 00 ,J 7 00 are attractors, 
the first one representing an accelerating solution. In the neighborhood of n = 1 the 
physics becomes more complex because V and Boo can acquire a stable character (node 
or focus) and C changes its stability. Finally, the phase space analysis shows that, 
for these type of theories, the sign of the curvature remains fixed once it has been 
set by the initial conditions. This is an important aspect of dynamics of the vacuum 
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Point 


n < 1/2 


1/2 < n < 1 


n > 1 


A 


repulsive 


repulsive 


repulsive 


Boo 


saddle 


attractive 


saddle 


c 

^oo 


saddle 


saddle 


repulsive 




attractive 


attractive 


attractive 


£ 

'-'OO 


saddle 


repulsive 


saddle 


T 

•> oo 


saddle 


saddle 


attractive 



model because, if the coordinate y of the point C is different from the sign of the initial 
condition for the curvature, the solution C is never reached. This means that C cannot 
be considered a global attractor for the vacuum case and an accelerating expansion 
phase is not guaranteed. 

4. The Matter Case 

When matter is present (as a perfect fluid), each of the cosmological equations gain a 
new source term containing the energy density. This modification has deep consequences 
because, as we can see from f|14H15p . the matter terms are coupled with a generic power 
of the curvature. Now, since in general the sign of the Ricci scalar is not fixed, these 
terms will not be defined for every real value of n. In particular we can allow only 
values of n belonging to the set of the relative numbers Z and the subset of the rational 
numbers Q which can be expressed as fractions with an odd denominator 

In this sense, the presence of matter induces a natural constraint through the field 
equations on i? n -gravity. In other words, if we suppose the field equations (|IJ) to be 
valid, not all the higher order theories of gravity are consistent in presence of standard 
matter. This difference between vacuum and non-vacuum physics is not present in GR 
but it is a common feature of higher order theories. 

The constraints on n also introduces technical problems because some of the results 
that we will present below require numerical calculations. This means that the result 
should be expressed in terms of the allowed set of values of n, but this is neither an easy 
nor a useful task. For this reason, we will work as if n was unconstrained, supposing 
that all the intervals that we devise are meant to represent the subset of allowed values 
of n within these intervals, and that their boundaries approximate the nearest allowed 
value of n. 

With this in mind, let us now derive the dynamics for the matter case. We again 
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start with (|15|) . and define the expansion normalized variables 

l-n), (45) 



6nH 2 

P 



3n X (n)H 2 R n - 1 
K " 



a 2 H 2 ' 

where z represents the contribution of standard matter in the form of the energy density 
weighted by a power of the Ricci curvature. In order to write down a dynamical equation 
for this variable, we need an expression for the time derivative of the energy density. Such 
an expression can be obtained from the Bianchi Identities T™ T ,u = which describe 
the conservation properties. As we have seen, if the gravitational field equations ((Zj) 
hold, it is possible to show 121] that by virtue of the Bianchi identities, whatever the 
form of the gravitational Lagrangian, the conservation equation for matter are simply 
given by (fTTjl . 

Differentiating (j45|) with respect to the logarithmic time, we obtain the autonomous 
system 



x' = x - 2x 2 + ( 2 + nx ^ y _ (2 + X )K + 2 - z(l + 3w) , (46) 
1 — n 

y' = ^ + 2v (^ + 2 + K ) , (47) 

n-1 \1 — n J 

z' = -3z(l + w) + 2z (2 + + Kj -xz, (48) 

K' = 2k(i + ^- + k) , (49) 



n 

with the constraint 

1 + x-y-z + K = . (50) 

As before, we can reduce the number of variables by substituting ()50|) into the other 
equations and recognizing that the spatial curvature equation is a linear combination of 
the other three. In this way, the final system becomes 

. 2(n — 2)y 9 xy 

x' = — — -2x- 2x 2 — + 1 + x)z - 3zw , 

n — 1 n — 1 

y = — ^-j- [(3 - 2n)x - 2y + 2(n - l)z + 2(n - 1)] , (51) 
z' = z ( 2z — 1 — 3x 3w 



n — 1 

Note that, as in the previous system, y = =^ y' = 0, but now also z = =^ z' = 
holds. Consequently the two planes y = and z = correspond to two invariant 
submanifolds. Since for z = (foTj) reduces to (|2*T|). we would be tempted to interpret 
the plane z = to be the vacuum invariant submanifold of the phase space. This is not 
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entirely correct. In order to understand this point, let us write the energy density p in 
terms of the variables (j45j). We have 

p = zy n - x H 2n . (52) 

Form this formula it is clear that when z = (and y ^ 0) the density is zero whatever 
the value of the parameter n. Instead if y — (and z ^ 0) the behavior depends on 
the value of n. In particular, the density is zero for y = and n > 1 and divergent for 
n < 1. If both x and z are equal to zero for n < 1 the behavior of p is undetermined 
and can be found only by direct substitution in the cosmological equations. 

As a consequence, for n > 1 both (x, y) and the (x, z) planes are invariant vacuum 
manifolds, but if n < 1 the vacuum submanifold is not necessarily compact because it 
is made by (x, y) and other (in principle) submanifolds of the (x, z) plane. 

The fact that z = and y = are invariant submanifolds also suggests that there 
are no orbits which can cross these planes: if we choose initial conditions yi = (z{ = 0), 
the equations (JoTj) imply that the orbits will never leave the plane y = (z = 0). As a 
consequence, if the evolution has initial conditions for which j/j ^ (zi ^ 0) they can 
approach the plane y — only asymptotically. Consequently, like in the vacuum case, 
no global attractor exists in these models because no orbit can, in general, cover all of 
the phase space. 



4-1. Finite Analysis 

Setting x' — 0, y' — 0, z' — and solving for x, y, z we obtain the fixed points: 

0, x -* 0, z -> 0}, K = -1, 

0, x -> -1, z -> 0}, K = 0, 

4n-5 2(n-2) ] , s 

x -> ^ -A 2^0 , K = 0, (53) 



4 = 


{y 


B = 


{y 


C = 


{* 


£> = 


{y 


£ = 


{y 


T = 


{y 


Q = 


{• 



2n - 1 ' 2n - 1 ' 

> 2(n - l) 2 , x -> -2(n - 1), z -> 0}, JsT = 2n 2 - 2n - 1 

0, x — > — 1 — 3w, z — > —1 — 3w}, if = —1, 

►0, x -> 1 -3w, ^ -> 2 -3w}, K = 0, 

(n-l)(4n- 12w-3) 3(n-l)(l + w) 



2n 2 n 
n(13 + 9tw) - 2n 2 (4 + 3w) - 3(1 + w) 



K = . 



2n 2 

Note that the fixed points A, B, C, V have the same coordinates (j22| found in the 
vacuum case, with the additional condition that the variable z has to be zero. 

Also in this case we have two different classes of fixed points, one independent of 
the value of the parameters n {A, 8,8,^) and one that changes its coordinates with n 

(c,v,g). 

There are also values of n for which two of the fixed points (C,Q) acquire an 
asymptotic character, namely 0, 1/2. However, n = 1/2 is not an allowed value in this 
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Table 5. Coordinates of the fixed points and solutions for i?"-gravity with matter. 



Point 



Coordinates (x, y, z) 



Scale Factor 



A 
B 
C 

V 

£ 
T 

Q 



[0,0,0] 
[-1,0,0] 

2(n-2) 4n-5 







2n-l ' 2n-l ' 

[2(l-n),2(n-l) 2 ,0] 

[-1 -3w,0, -l-3w] 
[1 - 3w, 0, 2 - 3iu] 

3(n-l)(l+w) (n-l)[4n-3(w+3)J 
n ' 2n 2 • 

n(13+9w)-2n 2 (4+3w)-3(l+«i) 
2n 2 



a = ao(t — to) 
oo(* - to) 1/2 (only for n = 3/2) 

(l-n)(2n-l) 

a = a t ™- 2 

a = a t if k = 
a = a (t — t ) 
ao(* - to) 1/2 (only for n = 3/2) 



a = ao t 3 '^™) 



Point 


Matter Density 




p = 




p = 


C 


p = 




p = 




p = for n > 1 or p divergent for n < 1 


T 


p = 


Q 


P = Pot- 2n , 




p = (-l)"3- n 2 2n ~ 1 n n (l + w)~ 2n x 




(An - 3(1 + w)) n - l [2n 2 (A + 3w) - n(13 + 9w) + 3(1 + w)) 



case and n = corresponds to a constant Lagrangian for the gravitational interaction. 
These values of n are not interesting for our purposes and are therefore excluded. 

Merging occurs for w in the Zel'dovic interval, excluding n = 1, for A and B at 
n = 5/4, for C and V at n = |(1 + v3) and n = |(1 — \/3), for C and when n is 
solution of the equation (8 + 6w)n 2 — (13 + 9w)n + 3(1 + u;) = 0. None of this values 
are allowed in the presence of matter, so we can consider these fixed points to be always 
distinct. 

The stability of the fixed points can be obtained, as usual, by linearizing the system 
(jf)T|) and are summarized in Tables (0 El and Because the phase space is three 
dimensional, the general dynamical behavior is not easy to visualize. We therefore do 
not give any graphical representation as in the vacuum case and focus only on the 
stability of the fixed points. 

The fixed point A has always a saddle-node character because the eigenvalues of 
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Point Eigenvalues 



A (-2,2,-1-3^) 
B (2,^,2 -3w) 

in— 5 2+An— 4n 2 — 3— 13n— 8n 2 — 3w+9nw— 6n 2 w 



c 



n-1 ' l-3n+2n 2 ' (ra-l)(2ra-l) 



V (-2 + n- v / 3n(3n - 4), -2 + n + v^^n - 4), 2n - 3(1 + w)) 
5 (-2,^=^,1 + 3^) 

^ (2,^^,3^-2) 

/ 3-2n+3w Pi(n,w)-y/P 2 (n,w) P 1 (n,w)+ y /P 2 (n,w) \ 
y In' 4n(n-l) ' 4n(n-l) J ' 

P^n, w) = (3(1 + w) + 3n((2n - 3)w - 1)) 

P 2 (n,iu) = (n-l)[4n 3 (8 + 3u0 2 -4n 2 (152 + 3w(55 + 18u0) 
+3n(l + w)(139 + 87w) - 81(1 + wf] 



Table 7. Stability of the fixed points for ^"-gravity with matter (dust or radiation). 
The term "spiral" has been used for pure attractive focus-nodes, the term "anti-spiral" 
for pure repulsive focus-nodes, and the term "saddle- focus" for unstable focus-nodes. 
The point Q has been treated separately because of the approximations used. 





n< |(1- v^) 




- y/3) < n < 


< n < 1/2 


1/2 < n< 1 


A 


saddle 




saddle 


saddle 


saddle 


B 


repellor 




repellor 


repellor 


repellor 


C 


attractor 




saddle 


saddle 


attractor 


V 


saddle 




attractor 


spiral 


spiral 


£ 


saddle 




saddle 


saddle 


saddle 


T 


saddle 




saddle 


saddle 


saddle 





1 < n < 5/4 


5/4 < n < 4/3 


4/3 <n< ±(1 + V3) 


n> ±(l + V3) 


A 


saddle 


saddle 


saddle 


saddle 


B 


saddle 


repellor 


repellor 


repellor 


C 


repellor 


saddle 


saddle 


attractor 


V 


spiral 


spiral 


attractor 


saddle 


£ 


saddle 


saddle 


saddle 


saddle 


T 


saddle 


saddle 


saddle 


saddle 



Cosmological dynamics of R n gravity 



18 



Table 8. Stability of the fixed points for i?"-gravity with stiff matter. The term 
"spiral" has been used for pure attractive focus-nodes, the term "anti-spiral" for pure 
repulsive focus-nodes, and the term "saddle-focus" for unstable focus-nodes. The point 
Q has been treated separately because of the approximations used. 



n<i(l- v / 3) |(l-v / 3)<n<0 0<n<l/2 1/2 < n< 1 1< n < £(11 + V^7) 



A 


saddle 


saddle 


saddle 


saddle 


saddle 


B 


saddle 


saddle 


saddle 


saddle 


saddle 


C 


attractor 


saddle 


saddle 


attractor 


repellor 


V 


saddle 


attractor 


spiral 


spiral 


spiral 


£ 


saddle 


saddle 


saddle 


saddle 


saddle 




saddle 


saddle 


saddle 


saddle 


saddle 





£(11 + y/Wf) < n < 4/3 


4/3 <n< KI + a/3) 


1(1 + ^3) <n<3/2 


n > 3/2 


A 


saddle 


saddle 


saddle 


saddle 


B 


saddle 


saddle 


saddle 


saddle 


C 


saddle 


saddle 


attractor 


attractor 


V 


spiral 


attractor 


saddle 


saddle 


£ 


saddle 


saddle 


saddle 


saddle 


T 


saddle 


saddle 


saddle 


repellor 



the linearized matrix do not depend on n and have an alternate sign whatever the value 
of w. 

The fixed point B is an unstable node for n < 1 and n > 5/4 if the standard matter 
is dust or radiation and a saddle otherwise. In the stiff matter case this point is always 
a saddle. 

The fixed point C behaves like a pure unstable node for l<n<5/4if the standard 
matter is dust or radiation and forl<n<£(ll + y/37) in case of stiff matter. A pure 
stable character is instead found for n < |(1 — \/3), 1/2 < n < 1 and n > |(l + ^3), 
whatever the value of the parameter w. For all the other values of n this point is a 
saddle-node. 

The fixed point V has eigenvalues that involve a square root of a function of n. 
This means that when this function is negative, the eigenvalues are complex and these 
points behave like focus-nodes. In particular, for n < and n > 4/3 this point is a 
pure attractive node if |(1 — v^3) < n ^ and 4/3^n<|(l + \/3) or a saddle-node 
otherwise. For < n < 4/3 the eigenvalues are complex and an analysis of the real 
parts show that the focus-nodes are always attractive. 

The fixed point £ is a saddle-node whatever the values of n and w. 

The fixed point T is a saddle-node for every value of n if the matter is dust or 
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Table 9. Stability of the fixed point Q. The term "spiral" has been used for pure 
attractive focus-nodes, the term "anti-spiral" for pure repulsive focus-nodes, and the 
term"saddle-focus" for unstable focus- nodes. 

77 < 0.33 0.33 < 77 < 0.35 0.35 < n < 0.37 0.37 < n < 0.71 0.71 < n < 1 

w = saddle saddle-focus saddle-focus saddle-focus saddle 
w = 1/3 saddle saddle saddle-focus saddle-focus saddle-focus 

w = 1 saddle saddle saddle saddle-focus saddle-focus 



l<n< 1.220 1.220 < n < 1.223 1.223 < n < 1.224 1.224 < n < 1.28 

w = saddle-focus saddle-focus saddle-focus saddle-focus 

w = 1/3 saddle-focus saddle-focus saddle-focus saddle-focus 

w = 1 saddle repellor saddle anti-spiral 



1.28 < ri < 1.32 1.32 < 77 < 1.47 1.47 < n < 1.50 n > 1.50 

w = saddle-focus saddle saddle saddle 

w = 1/3 saddle saddle saddle saddle 

w = 1 anti-spiral anti-spiral repellor saddle 



radiation. In the stiff matter case, this fixed point is an unstable node if n > |. 

The calculations of the eigenvalues of the fixed point Q cannot be performed in 
an exact way for a general n because they involve the solution of complete algebraic 
equations of order higher than two. Numerical calculations show that there are complex 
eigenvalues for 0.33 < n < 0.71 and 1 < n < 1.33 if w = 0, 0.35 < n < 1.28 if w = 1/3, 
0.37 < 77 < 1 and 1.224 < n < 1.47, if w — 1. For all these values of ra, behaves like a 
saddle focus. For the other values of n this point is always a saddle-node apart in case 
of stiff matter, for which 1.220 < n < 1.223 and 1.47 < 77 < 1.5. In this case we have a 
pure repulsive node % . 

Using the same procedure of the vacuum case, the coordinates of the fixed points 
(jo^|) can be used to find the behavior of the scale factor. The equation (J2l!j) generalizes 
to 

H = [xc + -zc-l^jH 2 . (54) 

If 77 1 and 

(n-l)(x c -z c -l)-y c ^0 (55) 

% Also in this case we will not investigate the behavior of this model in the bifurcations leaving this 
task for a future paper. 



Cosmological dynamics of R n gravity 20 
this equation can be integrated to give 

(yc \ 1 
1 — xc V zc \ ■ (56) 
n-1 J 

However, the characterization of the solution in this case requires also the derivation of 
the energy density p which is given by the (|52p. 

Since for z c = the solution (J5fi)l reduces to the the fixed points A, B, C, V are 
associated with the same solutions given in the vacuum case. The points C and V have 
y ^ and represent vacuum solution for every value of n and using the cosmological 
equations is it possible to show that this is also true for the solution associated with A 
and B. 

The equation (|56j) also allows one to find a solution for £. This point correspond 
to a vacuum Milne solution, but since the energy density in this point is divergent for 
n < 1 we consider this point physical only for n > 1. 

The solution for point JF, instead, cannot be found via the (jo^|) because the 
condition is violated. A direct resolution of (j^IJ) tells us that this point represents 
a solution only if p = and n = 3/2, but this value of the parameter is forbidden in 
this case. 

The fixed point Q is instead a new solution with respect to the vacuum case. The 
scale factor and the Hubble parameter, are in this case, 

2 11 

a = a (t -t ) 3(1+u,) (57) 
and the energy density is given by 

P = Pot~ 2n (58) 

with 

p = (-l) n x (n)3~ n 2 2n - 1 n n (l + w)~ 2n x (59) 
(An - 3(1 + w)) n - 1 \2n 2 (A + 3w) - n(13 + 9w) + 3(1 + w)} . (60) 

Thus Q represents a power law regime which, for n > 0, is an expanding solution with p 
decreasing in time, and for n < a contracting solution with p increasing in time. If n 
is close to 1, this solution can be thought of as an "almost- Friedman" solution. Instead, 
for n > |(1 + w) this solution corresponds to accelerating expansion. 

Unfortunately, since po is a function of n and w, p is not necessarily positive and 
consequently Q is not necessarily a physical point. The values of n for which Q is 
physical can be found by solving the inequality p > 0. We can distinguish two different 
cases: n belonging to the set of even integers or allowed rational numbers with an even 
numerator, which we call Af eV en, and n belonging to the odd integers or allowed rationals 
with odd numerator, which we call M dd- 

{n G Neven) in this case we require that 



(An - 3(1 + w))[2n 2 (A + 3w) - n(13 + 9w) + 3(1 + w)} > 
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This inequality is satisfied for 0.28 < n < 0.75 and for n > 1.35 if w = 0, for 
0.31 <n<landn> 1.29 if tu = 1/3, for 0.35 < n < 1.22 and n > 1.5 if w = 1. If 
n < 1 we should also add the constraint (4n — 3(1 + w)) 7^ to avoid divergences. 
However, the solution of this equation for physical choices of w have values of n 
which are not forbidden by the above inequalities, or do not correspond to the case 
n = 1 . As a consequence, the matter energy density always remains finite. 

(n G Modd) This case can be split in two different subcases 

n > for which 

[2n 2 (4 + 3w) - n(13 + 9w) + 3(1 + w)} < , 

which leads to 0.28 < n < 1.35 if w = 0, 0.31 < n < 1.29 if w = 1/3, 
0.35 < n < 1.22 if w = 1, the value n = 1 being excluded, 
n < for which 

[2n 2 (4 + 3w) - n(13 + 9w) + 3(1 + w)} > , 
leads to n < for w = 0, 1/3, 1. 



T , . , 13 + 9w ± V73 + 66w + 9w 2 . . . . 

it should be also pointed out that for n = ; r , the density is 

P 2(8 + 6w) ' 3 

zero. However these values are not allowed in the matter case, so the matter density is 
non-zero for all physical values of n. 

Using (|52jh the conditions above can be also read as conditions on the coordinates 
of G. In particular, the point G will represent a physical point only if 

zy > for ne Af even ; 

z < 0, n > r (61) 

„ for neModd- 
z > 0, n < 

Again, it is useful to define the deceleration parameter q and obtain the regions of 
the phase space in which correspond to accelerated evolution. We have 

= z-x — (62) 

n — 1 




y>0 1 (z-x)(n- 1) -y < 1 ' 

Thus, chosen n, the phase space is divided in two regions which correspond to 
accelerating and decelerating evolution. It is easy to check that the conditions (jfiHj) 
are never satisfied for the points A,B,V and for £,T (as we would expect from the 
solutions given above). For the point C the constraint on n are the same as (}33|) and 
for the point Q we have an accelerated evolution only if n > |(3 + 3a?) (expansion) and 
n < (contraction). 
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4-2. Asymptotic analysis 

As already discussed for the vacuum case, we complete the phase space analysis by 
investigating the asymptotic regime. The compactification of the phase space is achieved 
in the same way as the vacuum case, the difference being that now we now use spherical 
polar coordinates (r, 9, <fi) for which 

x — > r sm(9) cos(0), y — > r sm(9) sin(0), z — > r cos(6 ) ) , (64) 

where r G [0, oof, 9 G [0, 7r], cf> G [0,27r]. As in the vacuum case, the actual 
compactification is achieved transforming the radial coordinate r in 1Z. In the limit 
1Z — > 1 the dynamic equations (joTjl become 

ft' -> 1 {8(n - 1) cos 3 (#) - 2(n - 1) cos(#) sin 2 (fl)(cos(2#) - 3) - 4cos 2 (#) sin(#)x 
4(ri — 1J 

[3(n - 1) cos(0) + 2sin(0)] + sin 3 (#) [(9 - 8n) cos(</>) - cos(30) - 7sin(0) 

+ sin(30)]} (65) 

0/ _^ cos(0) sin(2g){2n sin(g) + 2cos(0)[(n - 1) cos(g) + sin(g)(sin(0) - cos(0))]} 
~* 4(n-l)(l-7i) ^ ^ 

^ 07 ^ (n _ 1) 1 (1 _^ ) { cos W sm(0)[(n - 1) cos(fl) + sin(^)(sin(0) - cos(0))]} . (67) 

This system has many interesting features. First of all, the above equations are 
parameterized by n but not by the barotropic factor w. This means that the asymptotic 
regime is independent of the nature of the perfect fluid that we choose. In addition, 
since there is no TZ dependence in the (|65p. the fixed points are determined only by 
the angular equations. Therefore, in order to find the asymptotic fixed points, we can 
neglect the 1Z' equation and deal only with the (9', <fi') system. 

Looking carefully at the system (|66H67|) we realize that it presents a non-trivial 
structure because the polynomials in the R.H.S are not prime to each other. This 
means that we have an infinite number of fixed points due to the zeros of the common 
factor. These points will define "fixed subspaces" of the total phase space in which every 
single point has, in principle, its own stability. 

Setting 9' = and = 0, we obtain the coordinates of the fixed points. Although 
we obtain 26 independent solutions, many of them are associated with the same fixed 
point on the sphere due to the periodic nature of the trigonometric functions in (j66H67j) 
and the peculiar topology of the sphere. The results are summarized in Table El 

We can distinguish a first set of eight ordinary fixed points whose coordinates do 
not depend on the value of n and are therefore common to all the R n models, two "fixed 
lines" LI and L2 again independent of n and other four fixed points on LI and L2 whose 
coordinates depend on n. 

The solutions associated with the first four (AxrL'oo) can be obtained using the 
same kind of procedure used in the vacuum case, we can obtain their solutions (see 
Table CU). 
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Table 10. Coordinates of the asymptotic fixed points and relative growing rates in 
non-vacuum i?"-gravity. Here "ind." represents an indeterminate form. 



Point 



{OA) 



Coordinates 



Growth Rates 

y/x z/x y/z 



A 

Boo 

c 

^00 

Voo 

'-'CO 

Goo 



(0,0) 

(7T,7r) 

(tt/2,0) 

(-7T/2,7r) 

(tt/4,0) 
(3tt/4,tt) 
(7r/2,7r/4) 
(tt/2,5tt/4) 



x,y 
x,y 
2/, 2 

X, 2 
X, z 

x,y 
x,y 



,z — > +00 
0,2^ —00 
, x — ► +00 



,x 



-00 



+00 , y — > 
-00 ,y -> 
+00 ,2^0 
—00 ,2^0 



00 

00 

1 ind. 

1 ind. 

1 

1 

1 00 
1 00 



Line LI 
Line L2 



(M/2) 
(Mtt/2) 



y, 2 — > +00 , rr — > 
y, z — > —00 , x — > 



00 00 tan(0) 
00 00 tan(0) 



J 1 

00 


(arctan(n - 


-1W2) 


V,z- 


* +00 




-> 


00 


00 


n — 


1 


J 2 


(— arctan(n 


-1),t/2) 


y,z- 


* +00 


,x - 


■» 


00 


00 


n — 


1 


J 3 


(arctan(n — 


1),3tt/2) 


y,z- 


■> —00 


,x - 


-> 


00 


00 


n — 


1 


J 4 


(— arctan(n - 


-1),3tt/2) 


y,z- 


■> —00 


,x- 


■* 


00 


00 


n — 


1 



Table 11. Coordinates, eigenvalues and value of r' of the ordinary asymptotic fixed 
points in non- vacuum i?"-gravity. 



Point 



{OA) 



Eigenvalues 9 — 



Behavior 



^4-00 
Boo 

F 

Goo 

7~Loo 



(0,0) 

(7T,7r) 

(tt/2,0) 

(-7T/2,7r) 

(tt/4,0) 
(3tt/4,tt) 
(tt/2,tt/4) 
(tt/2,5tt/4) 



[1,1] 
[-1,-1] 

2 ' 2(n-l) 

2 ' 2(l-n) 
y 7 ^ ny/2 



2(l-n) ' 2(n-l) 



2(n-l) ' 2(l-n) 



2 
-2 
-2 

2 

2 

VI 
2 

(l-2n)>/2 
2(l-n) 

(2ra-l)y / 2 
2(l-n) 





- A/'oo 


|AT 


- A/'oo 


|AT 


- A/'oo 




- ■A/'oo 




a = 




a = 



|A/"- A/^l = [Ci±Co|^|(*-*o)] 
|A/- A^oo I = [Ci±C |^|(*-to)] 



2n-l 
n-1 

2n-l 
n-1 
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Table 12. Stability of the ordinary asymptotic fixed points for i?"-gravity with matter. 
These results are independent of the value of w. 



Point 


n < 


< n < 1 


n > 1 


A 


repellor 


repellor 


repellor 


Boo 


attractor 


attractor 


attractor 


c 


saddle 


saddle 


saddle 


£>oc 


saddle 


saddle 


saddle 




saddle 


attractor 


saddle 


<*\ 30 


saddle 


repellor 


saddle 


Qoo 


saddle 


saddle 


saddle 


T~toc 


saddle 


saddle 


saddle 



For the points and J 7 ^, y — > and x, z — > ±oo at the same rate (x/z — > 1). For 
these points equation (J53j) reduces to if = which corresponds to the de Sitter solution 
a = a exp[A(£ — t )} where A is an integration constant. 

For the points Q^, Hoc, z — > and x, y — > ±oo at the same rate (x/y — > 1). In this 
limit the first and third equation of the system (|5T|) reduce to 

2 



x 



—x 



which admits the solution 



x 



Substituting this result in the (|53j) we obtain 

1 



n 



{t - 1 ) 



(68) 



(69) 



2n - 1 

as expected from the fact that for z = the system (j51|) reduces to (|21)1. 

The parameter n also determines the stability of the ordinary fixed points. As we 
can see from Table H2 all but the point and are unstable for every value of n. 
The point B^ is always an attractor and the point £qo is an attractor only if < n < 1. 

In addition to the ordinary fixed point, there is a second class of fixed points defined 
by two " fixed lines" LI with cf) = tt/2 and L2 with <fi = 3n/2 in the (9,<p) plane. The 
solution associated with the fixed lines can be obtained with the same method used 
above. We obtain 

2 



tn 



(70) 



For the points on LI and L2, one of the eigenvalues of the associated linearized matrix 
is null so that the system is structurally unstable. This means that the linearization 
is inappropriate and the local phase portrait is determined principally by the non- 
linear terms. In the presence of structural instability there is no simple way to devise 
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the stability properties of a fixed point and we are forced to study directly small 
perturbations around the fixed lines. Developing the (j66H67|) at first order about the 
generic point (vr/2, #0), we obtain 

89' 



n sin 



i ) cosltfoj 



n — 1 
(n — 1) cos(6> ) — sin( 



n 



This system admits the solutions 
Ainsm 2 ( 



86 



yp) cos^oj 

sin(^o) + (1 — n) cos(# ) 
(n — 1) cos(6> ) 



cxp 



n 



1) cosf 



AA 



n 



5i, 



Ai exp 



sm(0 o J 



1 



AT 



(71) 
(72) 

(73) 
(74) 



where the " ± " is connected with the sign of the perturbation 8<p and A% and B\ are 
integration constants. Now, if we look at the value of r' at the generic point (71/2,80) 
we obtain 

(n - 1) cos(#) - sin(0) 



f \4>=tt/2 



n 



(75) 



and discover that the stability of the fixed point depends only on the sign of the term 
(n — 1) cos(6 l o) + sin( " 



— . If M increases with time (i.e. the model is expanding) and 

is positive, the points on the fixed line (p = tt/2 are saddle nodes 

is negative, the points on 



n — 1 
the perturbation 

for every value of 9 and n. Conversely, if the perturbation 
the fixed line <fr = tt/2 are pure attractors if arctan(n — 1) < 9 < tt/2 for n > 1 and 
— 71/ 2 < 9 < arctan(n — 1) for n < 1 and are pure repellors otherwise. This behavior is 
reversed in case of contraction (M decreases with t) . 

Similarly, expanding the equations (|66H67|) at first order about the generic point 
(37r/2,#o) we obtain 

nsin 2 (6'o) cos(#o) 



89' = 6(f>- 



1 — n 
sin(6 l o) + (n 



I) cos( 



1 — n 

This system admits the solutions 

A 2 nsin 2 (# )cos(# ) 
89 = ± — — r — — exp 

sin(fc'o) + (n — 1) cos(#o) 

sin(6 l o) + (n — 1) cos 



sin(#o) + {n — 1) cos( 



71 



Bo 



8(f> = A 2 exp 



1 — 71 



(76) 
(77) 

(78) 
(79) 



where the " ± " is connected with the sign of the perturbation 5<ft and A2 and B2 are 
integration constants. Again, if we look at the value of r' at the generic point (vr/2, 9q) 

(n - 1) cos(0) + sin(6>) 



r \<j>=3ir/2 



(80) 



n — 1 

we discover that the stability of the fixed point depends only on the sign of the term 
(n — 1) cos(6' ) + sin(6> ) 



n—1 



If M increases with time (i.e. the model is expanding) the 
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perturbation 5(f) is positive, the points of the fixed line (f> = 3n/2 are saddle nodes for 
every value of 9 and n. Conversely, if the perturbation 5<f) is negative, the points of the 
fixed line (f> = 37r/2 are pure attractors if —it/2 < 9 < arctan(l — n) for n > 1 and 
arctan(l — n) < 9 < n/2 for n < 1 and are pure repellors otherwise. This behavior is 
reversed in case of contraction. 

Finally, there are four special points given by (±arctan(n — l),7r/2 + kir) with 
k = 0,1. These points belong to the lines LI and L2 and are associated with the 
solutions (|70|h The difference is that, in calculating their stability we find that both 
the eigenvalues are null. This situation is more complex than the previous one, and the 
stability analysis requires us to perturb the system at the second order. Unfortunately, 
the solution of this system turns out to be achievable only numerically and hence makes 
it impossible to give a general stability analysis. For this reason we will limit ourselves 
to giving only plots of the solutions for the point (tt/2, arctan(n — 1)) and for a set 
of specific values of n (see figures 9-14). It is clear that for all the values of n the 
perturbation on 9 approaches a constant. Such a behavior is due to the fact that the 
fixed lines are themselves invariant submanifolds and therefore small perturbations on 
9 induce the state to evolve on the line towards the nearest attractor. 

It is clear that the global behavior of the non-vacuum i?"-gravity cosmological 
models is not trivial. We have three different types of transient states which represent, 
Milne and "Friedman-like" solutions. As in the vacuum case, for values of n far from 
one the power law solution C is an attractor. When n approaches unity, other attractors 
i.e. T> come into play and the stability of C changes dramatically. Asymptotically two 
different sets of attractors are present. The first one is made up of ordinary points i.e. 
Boo and, for < n < 1, Soo] the second one by the arcs (</> = tt/2, < 9 < arctan(n — 1)) 
and {(f) = 37r/2, arctan(l — n) < 9 < 0) for n > 1 and arctan(n — 1) < 9 < and 
< 9 < arctan(l — n) for n < 1. 

The most interesting scenario is an orbit that approaches first the Q point and 
then the C solution. Clearly, the values of n should ensure Q is physical and not too 
different from the standard Friedmann solution in order be consistent with observational 
constraints. At the same time, n has to be chosen to ensure that C is associated with 
an accelerating solution in order to get a late acceleration phase. Moreover, since the 
planes y = and z = are invariant submanifolds and cannot be crossed by any orbit, 
we have to make sure that the y coordinate of the points C and Q have the same sign. 
Using the Tables given above, we can select the value of n for which this is realized: 



w — 1 



w = 0,1/3 




n e M £1 
n e Af c 



n E M e ,, 

neAf c 



even 



odd 



even 



odd 



(82) 



(81) 
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5. Discussion and Conclusions 

In this paper, we have given a detailed analysis of the behavior of i? n -gravity 
cosmological models using the dynamical system approach. Defining an effective fluid 
representing the higher order corrections to standard GR, it is possible to define a system 
of first order differential equations which are equivalent to the cosmological equations. 
In this way, using standard dynamical systems theory, we have obtained a qualitative 
description of the global and asymptotic behavior of these cosmologies. This allows 
to achieve a comprehensive (but qualitative) description of the global features of such 
models. 

Our investigation began with an analysis of the the pure gravity case in which the 
phase space is two dimensional. We found five fixed points, two of them (one of which 
is double) with coordinates independent of the parameter n and two whose position on 
the plane changes with n . These fixed points correspond to the four exact solutions 
listed in Table Q One of them, namely the solution associated with the point C, is 
particularly interesting because it can give rise to accelerated expansion. This solution 
has already been found by solving the field equations directly ^Qj and was compared 
with observational data. Our analysis shows that, for the values of n compatible with 
observations, the solution C is an attractor. 

However, the structure of the phase space shows that once we have set the initial 
condition, the sign of the Ricci scalar is preserved. In other words there is no way in 
which these models can evolve from a state R > to a state R < or vice versa. 
Such a transition would be impossible in GR because R is proportional to the energy 
density and a change in its sign would imply a violation of the energy conditions. In R n 
gravity, and in a general f(R) gravity, this is not necessarily true because of the more 
complicated relation between R and p. However, the main consequence of such a feature 
is that, in the vacuum phase space, no global attractor can exist and, in particular, C is 
an attractor only for specific initial conditions. 

The introduction of the matter term in the field equations induces strong constraints 
on the parameters. The value of n has to be either a relative or an odd denominator 
rational number to ensure self-consistency of the theory. This was expected: not all 
forms of non-minimal coupling are necessarily compatible with standard matter. Other 
constraints come from the fact that the matter energy density is a function of n and is 
not positive definite in general. 

The analysis of the dynamical system reveals that four of the seven fixed points 
found in this case are generalizations of the fixed point of the vacuum case and present 
the same stability properties. This implies that, for example the constraints on n found 
in ^E] hold also in this case. The last three fixed points are just transients that, 
providing w lies in the Zel'dovic interval, are physical only for specific values of n. The 
fixed point Q is particularly interesting and may be interpreted as corresponding to an 
"almost-Friedmann" transient phase in the evolution of the cosmological model. 

We found values of n for which the universe undergoes a transient almost- 
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Friedmann behavior, finally evolving into an accelerating phase. The existence of such a 
transient would not only explain the recently observed non constancy of the deceleration 
parameter but would also allow, in principle, a period of structure formation. Such 
a feature is relevant because it is not easy to explain how large scale structure is formed 
in accelerating universes. 

In the asymptotic regime i? n -gravity possesses number of fixed points which is the 
same or even bigger than the number of finite ones. In the vacuum case, we found 
that these fixed points represent essentially an exponential evolution. In two cases 
the character of the solution (expanding or contracting) depends on the value of the 
parameter n. For n in the interval suggested by the observations we have two attractors 
for n < 1 and only one for n > 1. 

In the matter case, the asymptotic structure is very rich. We have two different 
sets of critical points. One of them is made by single points and one is made by two 
continuous subsets of the phase space (the "fixed lines"). For the set of points (Ax,- 
(/oo) the asymptotic behaviors are of the same type of the vacuum case apart a pure de 
Sitter evolution when x and z diverge at the same rate (Soo-^oo)- For the fixed lines the 
behavior is of the same type of (Axr^oo)- 

The results presented in this paper have two key features. Firstly, the application 
of the dynamical system approach to higher order gravity cosmology seems to be a 
powerful tool, revealing global features of these models. Secondly, the results seem to 
indicate that the idea of curvature quintessence i.e. to model the accelerated expansion 
phenomenon as a high order gravitational effect, not only seems viable, but is also in 
good agreement with our present understanding of the history of the universe. 
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